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Abstract. The Maximum Likelihood Method is generalized to include effects 
important for UHECR applications. The new approach can incorporate source 
distance constraints implied by the observed CR energy and can allow for energy 
uncertainties, possible deflection in magnetic fields, multiple source types, and 
a spectrum of CR composition. It can be efficiently implemented and does 
not require the unphysical "isotropic" assumption for unidentified sources. The 
approach optimizes the utility of UHECR data to discriminate between source 
classes and can help constrain galactic and extragalactic magnetic fields. Aspects 
of the method are directly applicable in other contexts, such as TeV gamma ray 
astrophysics. 



1. Introduction 

The Maximum Likelihood (ML) method aUows for measured or modeled uncertainties 
in a dataset to be incorporated in an analysis. It allows data with different intrinsic 
resolution to be combined, and makes most efficient use of all available information. 
Most importantly, it avoids the disagreeable feature of binned analyses on small 
datasets: fluctuations in the conclusion depending on exactly where cuts are made 
and thus which events are included. The basic idea is to associate with each event 
in the dataset, a probability density for an event to be found with the observed 
characteristics (e.g., energy, arrival direction and possibly other properties such as 
composition indicators), for a given hypothesis as to how the observed UHECRs were 
produced and propagated to the detector: the type of sources, deflection properties 
of intervening magnetic fields, etc. The production hypothesis is characterized by 
some assumptions (e.g., that all UHECRs are protons - thus specifying the energy 
loss process - and that UHECR sources reside in galaxies found in certain catalogs) 
and a small set of parameters (e.g., the number of UHECRs which originate in each 
source catalog, the mean-square magnetic deflection per unit propagation distance, 
and perhaps a luminosity threshold for the source). The likelihood of finding the 
observed data-set is proportional to the product of the probability densities of the 
events making it up, and the most probable value of the parameters - under the given 
hypothesis - are those that maximize the likelihood. 

Up to now, applications of the maximum likehhood method to UHECR data[lll3| 
have considered only arrival direction correlations and angular resolution. Energy 
resolution was not included and the constraint on probable source distance from energy 
losses during propagation (referred to collectively below as the GZK constraint, for 
short) was implemented only crudely, by imposing a hard cut on source distance. 



Generalized Maximum Likelihood Method for Ultrahigh Energy Cosmic Rays 



2 



However these are crucial features of the data and the ML method needs to be 
generahzed to incorporate them. The framework to do so, and to allow for multiple 
source types, possible deflection in magnetic fields, and CR composition as well, is 
presented below. A computational strategy for efficient implementation is presented. 

The approach optimizes the utility of a given UHECR dataset (which can come 
from different experiments) to discriminate between source classes and can help 
constrain galactic and extragalactic magnetic fields. The method is presented in terms 
of UHECRs but it can also be used to good effect in simpler contexts, such as TeV 
gamma ray astrophysics. 

2. Maximum Likelihood Method 

Since UHECRs may be produced by several different mechanisms (e.g., the jets of 
powerful radio galaxies [H, Gamma Ray Bursts young magnetarsjsj, and Giant 
AGN FlaresQ), the analysis must allow for multiple source catalogs labeled k. The 
total number of catalogs being tested is iVcat and there are galaxies in the fcth 
catalog. Let Vk be the (a priori, unknown) number of CRs coming from objects in the 
kth catalog and vtot the total number of UHECRs in the dataset, 'Eki'k = vtot- 

Given the angular resolution for the ith event - deferring discussion of magnetic 
deflection until section 0] to simplify this introduction - one can write the probability 
density to find that event at any position in the sky if it were produced by some source 
candidate (labeled j). Call that probability density function Qi{xi — Sj,ai, ...), where 
Ui characterizes the angular resolution of the event, Xi is its arrival direction, and sj 
the direction of the source. Unless deflection is large, flat geometry is adequate and 
X, Sj, a are 2D vectors in the plane of the sky; for brevity, flat geometry is adopted and 
2D vector notation is suppressed except where essential. The ellipses in Q{x— sj ,cr, ...) 
denote information that is relevant when magnetic deflection is taken into account, 
deferred to section [4] below; Q has unit normalization 

j d^x Q{x,a, (1) 

For the case of no magnetic deflection and isotropic Gaussian resolution, 

Q(x.-.„.)^^exp(-i:^); (2) 

for this distribution, a cone of radius 1.5a contains 68% of the events. 

For each event, one needs the relative probability for it to come from a source at 
different redshifts. Due to GZK energy losses, this depends strongly on the observed 
energy E above about 50 EeV. It also depends on the assumed CR production 
spectrum at the source, since a harder spectrum and higher maximum energy mitigate 
the energy attenuation to some extent. A natural assumption, unless one has a large 
enough data-set to potentially have multiple events from individual sources or there is 
external information, is that every candidate source in a given catalog produces CRs 
with the same spectrum, apart from overall normalization. It is more convenient to 
work with e = In E, and we denote the probability density to observe a CR in the 
range e, e -\- de from a source at distance z + dz by F{e, z, X); A characterizes the 
energy resolution. For a fixed energy of the observed cosmic ray, F{e, Zj, X)/ F{e, Zji , A) 
is the relative likelihood that the CR comes from a source at Zj compared to one at 

Zj,. 
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For example, we could take the source spectrum to be a power-law, E~p, with 
Cmin < £ < Emax and p — 2 OT 2.5. Thc results of the analysis will depend on the choice 
of p and emax- The analysis should be done for several model spectra (here, values 
of p and Cmax) to estimate the systematic uncertainty in the ML method due to the 
spectrum. It may also be possible to treat p and emax as parameters to be determined 
in the analysis. As long as emin is low enough that it includes the minimum true 
arrival energy given the energy resolution and Sthresh, and low enough that energy 
losses in propagation do not move CRs from the most distant source out of the range 
e' > emin, the value of e„iin does not affect the result. The absolute normalization of 
F does not matter, but must be the same for each source. 

We will denote the true energy at Earth by E' and e' = In E' and define G'(e, z) to 
be the probability density for a CR reaching Earth from a source at distance z, z + dz 
to have a true In energy in the range e', e' + de' , given the assumed initial spectrum. 
If the energy uncertainty is just a log-normal energy resolution of width A, then 

F{e,z,\) = I de'G(e',z)-^exp (-^^^^) • (3) 

In order to appropriately call them probability densities, one can choose 

f de' G(e', z) = 1; f de F{e, z, A) = 1. (4) 

Note that one can prepare a lookup table for F{ei, Zj, A^) where e^, Aj are the In 
energy and its resolution, of the ith UHECR, and the Zj^s are the source redshifts. 
Then, the generalized ML (GML) analysis runs just as fast as it would ignoring energy 
resolution. If there is a known systematic overall energy shift for some set of events, 
those event energies would simply be shifted before doing the analysis. If it is desired 
to infer the most likely value of the energy shift of a particular experiment, then the 
energy of each event from that experiment would be multiplied by a factor, call it 
Kexpj and Kexp would be an additional parameter to be determined by maximizing 
the likelihood as discussed below; events from different experiments can be assigned 
different values of Kexp and these values determined independently if the data is 
adequate. 

With the above functions defined, the probability density for a given UHECR 

event to be found having thc observed parameters is just the sum of the probability 
density that it comes from each of the individual sources in the catalogs, which can 
be written as 

Viixi, e) = Sf^- (—) Sf4 wf^ Q{xi - Sj,ai, ...) F{ei, z„ K), (5) 

V ^tot / 

where w)j is the weight attributed to the jth source in the fcth catalog. Each catalog 
must have the same total weight, so that the relative weight of the catalogs is just 
I'kl^xot- This we enforce by the normalization condition 

EfiiWif = 1. (6) 

In the absence of any specific knowledge to the contrary, each source has equal a 
priori probability and thus equal weight. However it is clear that the probability of a 
particular source contributing to the dataset is proportional to the relative exposure, 
Tjj, of the UHECR dataset at position Sj. Furthermore, the flux from a source falls 
as z~'^, so in the absence of more information about individual sources, would 
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be taken proportional to r/j Zj ^. Additional information about the sources can be 

included in uij*^^ if available. For instance, sources in the IRAS catalog might be 
weighted with their infrared luminosity. Thus a typical choice for weights might be 



(fc) 



(k) 

where /j ' is the measured flux of the jth element of the kih catalog. Note that if 
the number of objects in a catalog is large, the normalization condition ^ implies 
that the weights associated with individual objects in that catalog are correspondingly 
reduced, as they should be. Multiple source candidates may fall within the angular 
domain of a given UHECR, in which case all of them would contribute according to 
their weights. 

Since the candidate source catalogs cannot in general be assumed to actually 
contain all the sources of the events, a "background" option must be included in 
the k catalogs in eqn Traditionally this has been taken to be an isotropic 

distribution, but on astrophysical grounds one expects the sources of UHECRs to 
be concentrated in regions of higher matter density. Thus a better treatment would 
use the 3D density distribution derived from a complete, unbiased catalog of galaxies. 
In practice obtaining such a density distribution is not trivial, due to incompleteness 
in the Galactic plane for optical surveys and the problem of too-low sampling for 
spatially more complete catalogs like HIPASS as used in ref. Q. 

The traditional method of allowing for "background" adds an isotropic term to 
eqn (IH) 

... (8) 

t'tot 

where R{x) is the relative acceptance in the direction x (normalized so that 
/ d'^x R{x) = 1). This approach distorts the results of the analysis and should only be 
used for qualitative or preliminary studies. For instance, in a study for which UHECRs 
are compared to known BL Lacs versus an isotropic source distribution, the analysis 
assuming isotropic background will return a larger number of BL Lac as sources than 
if the background were taken to be the galaxy distribution, even if the true UHECR 
sources are some other type of galaxy and not BL Lacs, simply because BL Lacs are 
clustered with galaxies. 

A better way to allow for the possibility that the chosen source catalog(s) do not 
contain all the sources of the UHECRs in the dataset, is to include a "background 
catalog" in the list of k possible source catalogs. The background catalog should 
be volume-limited and unbiased. The criterion for adequate sampling density for a 
background catalog is that some source in the catalog should be within an angular 
separation < ct of each UHECR, and within the redshift separation over which 
F{e, A, z) is slowly varying out to the largest contributing redshift, given the energy 
threshold of the UHECR dataset. If that condition cannot be satisfied with a given 
all-sky catalog, then the choices are to use a catalog such as 2MASS Redshift Survey 
(2MRS) and use its mask to apply a corresponding hard cut to the domain from 
which the CR dataset is taken, increase the energy threshold of the UHECR dataset 
to reduce the relevant range of z so it is low enough that the condition can be met, 
or a combination of the above. 

Note that it is in general incorrect to use an isotropic distribution to "complete" 
a catalog in either redshift or arrival direction. This applies to either the background 
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catalog or to a source-candidate catalog. As shown in J], the ML method 
automatically imposes the condition that the global distribution of matches to each 
catalog should represent the (weighted) global distribution of elements of each catalog. 
Thus a non-realistic global distribution for a catalog also distorts the i/^'s obtained. 
This is why it is for instance not correct in principle to fill in a "background" catalog 
like 2MRS with an isotropic distribution in the region of the mask, or complete a source 
catalog like the Veron-Cetty Veron one, with an isotropic distribution for distance 
> 100 Mpc. On the other hand, using an isotropic distribution as an approximate 
completion of a global under-sampling for the background catalog may be acceptable; 
such a procedure needs further study with simulations. 

3. Maximizing the Likelihood 

The next step is to introduce a likelihood measure for the entire dataset, and maximize 
with respect to the unknown parameters, to find the most likely values of these 
parameters. In the simplest case, these are simply the k — 1 independent i^k values. 
It is computationally desirable to work with the logarithm of the product of the Vi^s 
for each of the individual events: 

LM = InVi. (9) 

Assuming Gaussian angular resolution, each Q falls off very rapidly. For instance 
restricting to source candidates within a radius Scr of a cosmic ray captures all 
but 1% of the likely sources. Furthermore, source candidates that are left out 
have Q^O.OlQpeak compared to a perfectly ahgned source with the given angular 
resolution. With such a restriction on the angular region included for each cosmic ray, 
computation of the LM of a given dataset for fixed values of the unknown parameters 
is straightforward, and not very time-consuming even if the background catalog is 
large. Standard methods for multi-parameter maximization can be used to find the 
parameter values that maximize LM. The test for whether a given catalog is relevant 
or not, is to redo the calculation successively eliminating the catalog with the smallest 
Vk- If the LM changes by more than one unit, that catalog is a significant contributor 
to the signal and should be kept. The error ellipsoid in the multidimensional parameter 
space can be determined by finding the parameter values that give LM values one, 
two, ... units lower than the maximum. Note that the parameters Vk are not in general 
integers - they are the mean in the sense of a Poisson distribution of the number of 
events contributed by fcth source class. 

The analysis can be extended to larger number of unknown parameters by 
first determining the essential catalogs and values of key parameters. Then, more 
parameters can be added, maximizing LM with respect to all parameters but now 
keeping the domain of variation of the original parameters small. 

4. Magnetic deflections 

Charged UHECRs are deflected during their propagation between source and Earth, 
but extragalactic and Galactic magnetic fields (EGMF and GMF) are poorly 
constrained. It may be possible to use UHECR correlations to constrain these fields 
and at the same time improve the source analysis discussed above. That is where 
the ellipses in Q{x — Sj,a, ...) come in. When magnetic deflections are included, Q 
depends not only on a; — sj and fi, but also on E, the absolute arrival direction x, and 
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additional parameters characterizing the magnetic deflections. To what extent it wiU 
be possible to both identify the sources and constrain the magnetic deflections depends 
on how rare and distinctive UHECR sources are and whether there are multiple source 
types. If sources are sparse and clearly identifiable, e.g., powerful AGNs, then multiple 
events from single sources will both strengthen the identiflcation of their source and 
constrain the magnetic deflection and dispersion along that line-of-sight. However if 
virtually every galaxy can produce UHECRs (e.g., at the birth of a magnetar) but 
such events are rare and the mean number of events from a single source is <C 1, it 
will be diflicult to use UHECRs to map the GMF. 

In the limit of small deflections, the effect of a large scale magnetic field can be 
described in terms of a displacement function, giving the angular displacement 

due to the Lorentz force per unit CR rigidity (energy/Z) experienced as a CR traverses 
the field. For large scale fields and high enough CR energy, the trajectory in the 
Galaxy can be approximated as almost linear and f3 varies slowly with arrival direction. 
Measuring UHECR energies in EeV, distances in kpc, angular separations in degrees, 
and magnetic fields in fiG, the deflection of a CR traversing a distance L through the 
GMF along the line of sight I is [3 Z E where 

/3 = 630- / fx Bdl. (10) 



L . 

Under these assumptions, the coherent displacement of events from a single source 
can be described by a common /3. The effect of magnetic deflection is incorporated in 
the GML analysis by changing Xi to Xi — p{xi) Zi Ei, so that Q becomes a function of 
energy. In practice, one would adopt some simple model of the GMF and maximize 
the LM with respect to its parameters. If the impact of energy resolution on arrival 
direction correlations is to be fully taken into account, the integral over e must include 
Q as well as G. 

In addition to coherent deflections by the GMF, propagation through turbulent 
random fields may be an important effect. In the limit of sufficiently many small, 
randomly oriented deflections, the smearing of a cosmic ray from a pointlikc continuous 
source due to random fields is described by a 2-d Gaussian probability distribution of 
width 

ol[E^D)^[^-) ^'-^^ , (11) 

where D is the distance to the source and A is the characteristic length scale of the 
turbulent fields (expected to be ~ Mpc for extragalactic but < 100 pc for Galactic 
random fields). The expression pT|) is easily understood in the case of randomly 
oriented magnetic domains of size A and strength _B, where the magnitude of the 
defiection with respect to the direction of motion in each domain is 89 = X/Rl 
and i?i = 1.08 kpc ^g'^"^^ . The total number of deflections is = D/X and 

the mean total deflection obeys A6'^ — N66^/3. The factor 1/3 compared to the 
familiar random walk formula was derived in [8| for the case of x-ray scattering off 
dust, and arises because A9 is the angle between the UHECR velocity vector and 
a vector pointing to the source, both of which change as the UHECR propagates. 
Isotropy implies < >— ^ < B^ >. The general expression for {B^X) as a 
weighted mean of B^X over the trajectory is given in Equation pT|) is applicable 
if a//V ^ \/D/ < A > ^1. If multiple events come from the same source, their 
trajectories must sample different magnetic domains for them to have independent 
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deflections, requiring that the difference in their displacements must be larger than a 
characteristic domain size, i.e.. 



^>(A), (12) 



Ei Eii 

for all events i, z' in the cluster under consideration. 

Since the experimental angular resolution of the ith event is approximately a 
symmetric 2-d Gaussian (of width (7^), the net result of Gaussian magnetic smearing 
is to increase the effective resolution of that event in computing Q for the jth source, 
to (Ji^cff — v^crf + (7B{Ei, Zi, DjY. For the simplest treatment in a correlation analysis 
between UHECRs and catalogs as discussed above, one could take 



with /3o a single universal parameter to be determined by maximizing the LM . Note 
that {B'^X) need not be constant for the description (fT^ to be adequate, as long 
as the distribution of random magnetic deflections is Gaussian; the inferred value of 
/3o would be the mean of the distribution. A more sophisticated treatment could 
be to assume the value of (B^XD) for the jth source candidate is proportional 
to some measure of the energy density of intervening IGM, e.g., the integrated 
luminosity of background galaxies along the line-of-sight to j - this would correspond 
to the reasonable assumption that the energy density in random magnetic fields is 
proportional the the luminosity of the galaxy. In this case, the proportionality constant 
could be determined by maximizing LM. A virtue of the GML method is that as the 
data-set increases, more complex hypotheses for the EGMF can be used. 



5. Composition 

In order to carry out the magnetic deflection analysis one needs to know or make an 
assumption about the charges of the UHECRs. As argued by [l3| it is likely that the 
highest energy cosmic rays have a pure or bimodal composition, due to the fragility of 
light and intermediate mass nuclei to photodisintegration during propagation. Thus 
it is natural to assume that all events are protons or Fe. However as UHECR 
reconstruction techniques improve, it will become possible to give individual UHECR 
events a probability distribution over primary particle type and hence charge, based 
on measured properties of the individual CR shower. If this information is available, 
eqn ^ is replaced by 

V. (x, , 6) ^ E£t ^ Sf^^ K,, [x, , e) (14) 

^tot 

where now 

Pi,j{xi,ei) ^J:z Q ' Q{xi- Sj,ai,Z,...)Fz{ei,Zj,Xi), (15) 
where Q is the probabihty for the zth event to have charge Z, normahzed to 

CP' = 1- (16) 

The energy loss during propagation depends on the charge of the UHECR, so F is 
generalized to Fz- Similarly, in the presence of magnetic deflection, Q depends on Z 
via eqns (jlOllip . 
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Prior to the time that reconstruction methods have improved to the point that 
individual events can be assigned Q values reliably, one can use eqn to constrain 
the fraction of particles of charge Z in the dataset. For instance if the highest energy 
cosmic rays are assumed to be a mixture of protons and iron, one would take (^^^"^ — fp, 

(26) ( Z) 

Q ='^-fp, aU other Q = 0- and maximize the LM with respect to fp. Clearly, 
the success of such an effort depends on having a large high energy dataset and source 
catalogs that contain the sources of most of the events. 



6. Application to clusters of events 

The GML method developed here works whether or not there are multiple events from 
a single source. However when a possible cluster is identified, as for instance the 4-5 



events of the "Ursa Major" cluster in the combined AGASA-HiRes datasetpjj, it is 
no longer necessary to correlate the members of the cluster to a candidate source, 
allowing for a self-contained analysis. Rather than using eqn ([5]) and maximizing the 
LM with respect to the i/^'s for the different source catalogs, one introduces a source 
position Xs, redshift Zg, and number Vs of events from that source. Like the Vk's in 
eqn ([5]), is not in general an integer. One maximizes the LM to determine the most 
likely source position and redshift, and the most likely number of events coming from 
the source. The alternative to coming from the source now being that they come from 
the background which as discussed above might be approximated by a catalog, b: 

T^ii,Xi^ 6j, Xs: Zs^ — I | Q{Xi Xg^ Cfi^ ...) -E{^i: Zg^ A^) 

V ^tot / 

1 - Efii wf^ Q{x, - s,,a^, ...) F{e,,z„{m 



Vtot 



The LM is formed in principle as in eqn ([9]), but the computational speed is 
greatly increased by separating the computation of the sum over UHECR events in 
eqn ([9|) into two parts, those which are near and those which are far from the cluster. 
The InViS of UHECRs which are too far in angle to be in the cluster, have a negligible 
Q value in the first term (proportional to Vg) in eqn (|17p . Furthermore, the factor 
multiplying (1 — ■^^) in the second term of (|17[) (the sum on j) is independent of Xg 
and Vg, the parameters to be varied. Thus the contribution of distant UHECRs to 
Yii InVi amounts to a piece that can be computed once and for all, plus Vi^x ln{l — -^^), 
where Vf^i is the number of distant events. Thus the maximization of LM over Xg and 
Ug is completely tractable computationally. 

The magnetic deflection and dispersion of the cluster events can be described 
by 3 parameters, (3 and E* , for a total of 7 parameters including Xg, Zg and Vg. A 
simpler version of this analysis for the Ursa Major cluster, without the energy-distance 
constraint and using the traditional isotropic background, was reported in [ll|, IT^ . 
The significance of the cluster is similar for either 4 or 5 events: ~ 0.2% probablity 
of being a chance fluctuation. Even with only 4-5 candidate events in the cluster, 
the results for /3 and E* are robust to variations in other parameters. The coherent 
magnetic deflection is negligible and {B'^\)D « 8nG^ Mpc^ (assuming Z = 1) (llj . 

The treatment of F described in section [2] above, using a power-law source 
distribution, applies only for a continuous source. A point emphasized in [ll| is that 
if a cluster of events comes from a bursting source, whose duration is short compared 
to the typical delay time from magnetic deflections, then the observed spectrum is 
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peaked rather than power-law and this must be taken into consideration. An analysis 
of the Ursa Major cluster in the spirit of the method outlined here, but applicable 
also for a bursting source, will be presented elsewhere. 

7. Discussion and Conclusions 

We have presented a general and practical way to maximize the utility of the 
ML approach for UHECR correlation studies. Its greatest immediate value is in 
incorporating the source distance information following from the GZK energy loss 
phenomenon. In first applications a simple composition could be assumed, but the 
method is general enough to allow for a distribution of possible charges of the UHECR 
and eventually a probability distribution for the charge of individual events when 
event-by-event composition indicators are improved. The method can be applied 
directly to the simpler case of TeV gamma ray data, to include the effects of energy 
loss from production in interactions with diffuse background radiation. A rescaling 
factor for the level of the background radiation in wavebands where it is poorly 
measured could be taken as a parameter to be determined by maximizing the LM. 

At present, the published high quality data on UHECRs consists of 57 events 
from AGASA with energies above 40 EeV |l3j . 271 events from HiRes above 10 
EeV (energy-ordering but not individual energies are available, G. Thompson private 
communication), and 27 events from Auger with energies above 56 EeV[l3|- (The 
energy and angular resolution of earlier experiments is lower and more uncertain so 
they have only marginal utility in a ML analysis.) With ~ 400 CRs whose typical 
resolution element is < 10 sq-degrees and a total aperture containing ~ 40k sq-degrees, 
and future analyses most likely oriented to the higher energy range to get the benefit 
of the GZK horizon, the dataset is small enough that the probability of two UHECRs 
from an isotropic dataset falling by chance within one characteristic angular resolution 
element is fairly small and the methods presented here should be quite powerful. 
Simulations to test its effectiveness under different scenarios of source density and 
catalog completeness are needed. 

The Auger collaboration has reported a statistically signficant correlation between 
the highest energy UHECRs and Veron-Cetty Veron galaxies with z < 0.018 within 



3.2°[15l Il4l|. An initial ML study applying the generalized method proposed here to 
the 27 published events, might assume a pure proton composition and allow for an 
overall energy rescaling parameter and a single magnetic smearing parameter as in 
eqn (jlip . This should have more discriminating power in a search for the sources 
of UHECRs, as a result of employing a more realistic background catalog and not 
requiring a hard cut in angular separation and Zmax as in the original, binned Auger 
analysis. Applied to the full Auger dataset including events with lower energy, one 
would expect that the number of events Vk attributed to the fcth candidate source 
catalog, and the significance of the correlation with it, should stay approximately 
constant as the energy threshold is lowered, with the number of events attributed to 
the "background" catalog increasing as the dataset does. This is because the "utility" 
of the added lower energy events to discriminate between source and background 
catalogs drops, as the number of candidate or background sources within a few-cr cone 
about the UHECR direction becomes large with increasing GZK horizon. Thus the 
method should be insensitive to the energy threshold of the UHECR dataset, as long 
as the redshift-completeness of the source and background catalogs are comparable; 
otherwise, artifacts can be introduced. Simulations to confirm this are underway. 
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